## To run the code, you need to execute the "replication_main_results" file.

###1. Bayesian factorial analysis: scores by dimensions----- 

##D1: Autonomy (Figure A1)

ggs_caterpillar(sample.post.agencies1, line = 0) +
  scale_x_continuous(name = "Autonomy",
                     breaks = seq(-4, 4, 1),
                     limits = c(-4, 4)) +
  scale_y_discrete(name = "Agencies") +#, #labels = levels(sample.post.agencies1$code))+ 
  theme_linedraw() +
  theme(axis.title.x = element_text(size = 15, face = "italic"),
        axis.title.y = element_text(size = 15, face = "italic"),
        panel.background = element_blank(),
        panel.grid.major = element_line(colour = "grey50", size = 0.1),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.grid.major.x = element_blank())

ggsave("gg_autonomy.png", width = 9, height = 6, dpi = 300)

##D2: Accountability (Figure A2)

ggs_caterpillar(sample.post.agencies2, line = 0) +
  scale_x_continuous(name = "Accountability",
                     breaks = seq(-3, 3, 1),
                     limits = c(-3.25, 3.25)) +
  scale_y_discrete(name = "Agencies") +
  theme_linedraw() +
  theme(axis.title.x = element_text(size = 15, face = "italic"),
        axis.title.y = element_text(size = 15, face = "italic"),
        panel.background = element_blank(),
        panel.grid.major = element_line(colour = "grey50", size = 0.1),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.grid.major.x = element_blank())

ggsave("gg_accoun.png", width = 9, height = 6, dpi = 300)

##D3: Auth.'s enforcement (Figure A3)

ggs_caterpillar(sample.post.agencies3, line = 0) +
  scale_x_continuous(name = "Authorities' enforcement",
                     breaks = seq(-3, 3, 1),
                     limits = c(-3.25, 3.25)) +
  scale_y_discrete(name = "Agencies") +
  theme_linedraw() +
  theme(axis.title.x = element_text(size = 15, face = "italic"),
        axis.title.y = element_text(size = 15, face = "italic"),
        panel.background = element_blank(),
        panel.grid.major = element_line(colour = "grey50", size = 0.1),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.grid.major.x = element_blank())

ggsave("gg_enforce.png", width = 9, height = 6, dpi = 300)

##D4: Consumers awareness (Figure A4)

ggs_caterpillar(sample.post.agencies4, line = 0) +
  scale_x_continuous(name = "Consumers Awareness",
                     breaks = seq(-3, 3, 1),
                     limits = c(-3.25, 3.25)) +
  scale_y_discrete(name = "Agencies") +
  theme_linedraw() +
  #ggtitle("Figure 5. General agencies' consumer awareness") +
  theme(axis.title.x = element_text(size = 15, face = "italic"),
        axis.title.y = element_text(size = 15, face = "italic"),
        panel.background = element_blank(),
        panel.grid.major = element_line(colour = "grey50", size = 0.1),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black"),
        panel.grid.major.x = element_blank())

ggsave("gg_conawe.png", width = 9, height = 6, dpi = 300)


###2. Correlation matrix (Table A3)----
library(openxlsx)

attach(dfgen_scatter)

CORdfgen_scatter <- data.frame(Mean.D1,Mean.D2,Mean.D3,Mean.D4)
cor_dfgen_scatter <- cor(CORdfgen_scatter, use = "pairwise.complete.obs")

write.xlsx(cor_dfgen_scatter, "CORRDIM_GENERAL.xlsx")

###3. Export dendogram plot (Figure A5)----

library("dendextend")

attach(dfgen_scatter)

#Plot
plot(df.dend, cex = 0.4)

png("dendrogram.png", 
    width = 1250, height = 1300)  
plot(df.dend, cex = 0.4)  
dev.off()

###4. Perform k-means clusters analysis with 4 groups (Figure A6)----
df_dist2 <- df_dend %>%
  dplyr::select(Mean.D1, Mean.D2, Mean.D3, Mean.D4) 

###Iterate the following code until having cluster 1 = 10, cluster 2 = 6,
###cluster 3 = 8, cluster 4 = 14

k <- kmeans(df_dist2, centers = 4)
table(k$cluster)

df_dend$k.clust <- k$cluster
df_dend$k.clust <- as.factor(df_dend$k.clust)

identical(rownames(df_dend), rownames(df_dist2))

k_mean <- df_dend %>%
  group_by(k.clust) %>%
  summarise(
    Autonomy = mean(Mean.D1, na.rm = TRUE),  
    Accountability = mean(Mean.D2, na.rm = TRUE),  
    Enforcement = mean(Mean.D3, na.rm = TRUE),  
    C.Awareness = mean(Mean.D4, na.rm = TRUE),
    .groups = "drop"
  )

k_long <- gather(k_mean, key = "Variable", 
                 value = "Value", -k.clust)

table(k_long$Variable)
k_long$Variable <- dplyr::recode_factor(k_long$Variable, 
                                        "Autonomy" = "Autonomy",
                                        "Accountability" = "Accountability",
                                        "Enforcement" = "Enforcement",
                                        "C.Awareness" = "Consumers\nawareness")

table(k_long$k.clust)
k_long$k.clust <- dplyr::recode_factor(k_long$k.clust,
                                       "2" = "Cluster 1",
                                       "1" = "Cluster 2",
                                       "3" = "Cluster 3",
                                       "4" = "Cluster 4")

ggplot(k_long, aes(x = Variable, y = Value, fill = "Group")) +
  geom_bar(stat = "identity", position = "dodge", alpha = 0.7, color = "#FFB6C1") +
  labs(x = "Dimensions", y = "Means of clusters' dimensions") +
  facet_wrap(~k.clust, ncol = 1) +
  theme_minimal() + 
  ylim(c(-2, 2)) +
  theme(legend.position = "none",
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        axis.text.x  = element_text(size = 13),
        axis.text.y  = element_text(size = 13),
        strip.text = element_text(size = 14))

ggsave("clust_kmeans.png", width = 13, height = 8, dpi = 300)


###5. Regressions' table----

summary(mc1)

modelsummary(list("Cluster 1" = mc1,
                  "Cluster 2" = mc2,
                  "Cluster 3" = mc3,
                  "Cluster 4" = mc4), 
             stars = TRUE,
             estimate  = c("{estimate} ({std.error}){stars}"),
             statistic = NULL,
             add_rows = NULL,
             coef_rename = c("(Intercept)" = "Intercept",
                             "inst_designContinental" = "Continental",
                             "inst_designEastern Europe" = "Eastern Europe",
                             "inst_designNordic" = "Nordic",
                             "inst_designSouthern Europe" = "Southern Europe",
                             "gdp_growth" = "GDP Growth (%)",
                             "gdp_percap" = "GDP per cap.",
                             "trade_gdp" = "Trade (GDP %)",
                             "tax_rev" = "Tax Revenue (GDP %)",
                             "inflation_gdp" = "Inflation (%)",
                             "euro" = "Euro (=1)",
                             "eu_membership" = "Years of EU membership",
                             "regime" = "Type of Regime (1-10)"),
             output = "Models1+Controls.docx")


